home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fit / linear.c next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  7.9 KB  |  347 lines
  1. /* fit/linear.c
  2.  * 
  3.  * Copyright (C) 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <gsl/gsl_errno.h>
  22. #include <gsl/gsl_fit.h>
  23.  
  24. /* Fit the data (x_i, y_i) to the linear relationship 
  25.  
  26.    Y = c0 + c1 x
  27.  
  28.    returning, 
  29.  
  30.    c0, c1  --  coefficients
  31.    cov00, cov01, cov11  --  variance-covariance matrix of c0 and c1,
  32.    sumsq   --   sum of squares of residuals 
  33.  
  34.    This fit can be used in the case where the errors for the data are
  35.    uknown, but assumed equal for all points. The resulting
  36.    variance-covariance matrix estimates the error in the coefficients
  37.    from the observed variance of the points around the best fit line.
  38. */
  39.  
  40. int
  41. gsl_fit_linear (const double *x, const size_t xstride,
  42.         const double *y, const size_t ystride,
  43.         const size_t n,
  44.         double *c0, double *c1,
  45.         double *cov_00, double *cov_01, double *cov_11, double *sumsq)
  46. {
  47.   double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;
  48.  
  49.   size_t i;
  50.  
  51.   for (i = 0; i < n; i++)
  52.     {
  53.       m_x += (x[i * xstride] - m_x) / (i + 1.0);
  54.       m_y += (y[i * ystride] - m_y) / (i + 1.0);
  55.     }
  56.  
  57.   for (i = 0; i < n; i++)
  58.     {
  59.       const double dx = x[i * xstride] - m_x;
  60.       const double dy = y[i * ystride] - m_y;
  61.  
  62.       m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
  63.       m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
  64.     }
  65.  
  66.   /* In terms of y = a + b x */
  67.  
  68.   {
  69.     double s2 = 0, d2 = 0;
  70.     double b = m_dxdy / m_dx2;
  71.     double a = m_y - m_x * b;
  72.  
  73.     *c0 = a;
  74.     *c1 = b;
  75.  
  76.     /* Compute chi^2 = \sum (y_i - (a + b * x_i))^2 */
  77.  
  78.     for (i = 0; i < n; i++)
  79.       {
  80.     const double dx = x[i * xstride] - m_x;
  81.     const double dy = y[i * ystride] - m_y;
  82.     const double d = dy - b * dx;
  83.     d2 += d * d;
  84.       }
  85.  
  86.     s2 = d2 / (n - 2.0);    /* chisq per degree of freedom */
  87.  
  88.     *cov_00 = s2 * (1.0 / n) * (1 + m_x * m_x / m_dx2);
  89.     *cov_11 = s2 * 1.0 / (n * m_dx2);
  90.  
  91.     *cov_01 = s2 * (-m_x) / (n * m_dx2);
  92.  
  93.     *sumsq = d2;
  94.   }
  95.  
  96.   return GSL_SUCCESS;
  97. }
  98.  
  99.  
  100. /* Fit the weighted data (x_i, w_i, y_i) to the linear relationship 
  101.  
  102.    Y = c0 + c1 x
  103.  
  104.    returning, 
  105.  
  106.    c0, c1  --  coefficients
  107.    s0, s1  --  the standard deviations of c0 and c1,
  108.    r       --  the correlation coefficient between c0 and c1,
  109.    chisq   --  weighted sum of squares of residuals */
  110.  
  111. int
  112. gsl_fit_wlinear (const double *x, const size_t xstride,
  113.          const double *w, const size_t wstride,
  114.          const double *y, const size_t ystride,
  115.          const size_t n,
  116.          double *c0, double *c1,
  117.          double *cov_00, double *cov_01, double *cov_11,
  118.          double *chisq)
  119. {
  120.  
  121.   /* compute the weighted means and weighted deviations from the means */
  122.  
  123.   /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */
  124.  
  125.   double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;
  126.  
  127.   size_t i;
  128.  
  129.   for (i = 0; i < n; i++)
  130.     {
  131.       const double wi = w[i * wstride];
  132.  
  133.       if (wi > 0)
  134.     {
  135.       W += wi;
  136.       wm_x += (x[i * xstride] - wm_x) * (wi / W);
  137.       wm_y += (y[i * ystride] - wm_y) * (wi / W);
  138.     }
  139.     }
  140.  
  141.   W = 0;            /* reset the total weight */
  142.  
  143.   for (i = 0; i < n; i++)
  144.     {
  145.       const double wi = w[i * wstride];
  146.  
  147.       if (wi > 0)
  148.     {
  149.       const double dx = x[i * xstride] - wm_x;
  150.       const double dy = y[i * ystride] - wm_y;
  151.  
  152.       W += wi;
  153.       wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
  154.       wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
  155.     }
  156.     }
  157.  
  158.   /* In terms of y = a + b x */
  159.  
  160.   {
  161.     double d2 = 0;
  162.     double b = wm_dxdy / wm_dx2;
  163.     double a = wm_y - wm_x * b;
  164.  
  165.     *c0 = a;
  166.     *c1 = b;
  167.  
  168.     *cov_00 = (1 / W) * (1 + wm_x * wm_x / wm_dx2);
  169.     *cov_11 = 1 / (W * wm_dx2);
  170.  
  171.     *cov_01 = -wm_x / (W * wm_dx2);
  172.  
  173.     /* Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2 */
  174.  
  175.     for (i = 0; i < n; i++)
  176.       {
  177.     const double wi = w[i * wstride];
  178.  
  179.     if (wi > 0)
  180.       {
  181.         const double dx = x[i * xstride] - wm_x;
  182.         const double dy = y[i * ystride] - wm_y;
  183.         const double d = dy - b * dx;
  184.         d2 += wi * d * d;
  185.       }
  186.       }
  187.  
  188.     *chisq = d2;
  189.   }
  190.  
  191.   return GSL_SUCCESS;
  192. }
  193.  
  194.  
  195.  
  196. int
  197. gsl_fit_linear_est (const double x,
  198.             const double c0, const double c1,
  199.             const double c00, const double c01, const double c11,
  200.             double *y, double *y_err)
  201. {
  202.   *y = c0 + c1 * x;
  203.   *y_err = sqrt (c00 + x * (2 * c01 + c11 * x));
  204.   return GSL_SUCCESS;
  205. }
  206.  
  207.  
  208. int
  209. gsl_fit_mul (const double *x, const size_t xstride,
  210.          const double *y, const size_t ystride,
  211.          const size_t n, 
  212.              double *c1, double *cov_11, double *sumsq)
  213. {
  214.   double m_x = 0, m_y = 0, m_dx2 = 0, m_dxdy = 0;
  215.  
  216.   size_t i;
  217.  
  218.   for (i = 0; i < n; i++)
  219.     {
  220.       m_x += (x[i * xstride] - m_x) / (i + 1.0);
  221.       m_y += (y[i * ystride] - m_y) / (i + 1.0);
  222.     }
  223.  
  224.   for (i = 0; i < n; i++)
  225.     {
  226.       const double dx = x[i * xstride] - m_x;
  227.       const double dy = y[i * ystride] - m_y;
  228.  
  229.       m_dx2 += (dx * dx - m_dx2) / (i + 1.0);
  230.       m_dxdy += (dx * dy - m_dxdy) / (i + 1.0);
  231.     }
  232.  
  233.   /* In terms of y =  b x */
  234.  
  235.   {
  236.     double s2 = 0, d2 = 0;
  237.     double b = (m_x * m_y + m_dxdy) / (m_x * m_x + m_dx2);
  238.  
  239.     *c1 = b;
  240.  
  241.     /* Compute chi^2 = \sum (y_i -  b * x_i)^2 */
  242.  
  243.     for (i = 0; i < n; i++)
  244.       {
  245.     const double dx = x[i * xstride] - m_x;
  246.     const double dy = y[i * ystride] - m_y;
  247.     const double d = (m_y - b * m_x) + dy - b * dx;
  248.     d2 += d * d;
  249.       }
  250.  
  251.     s2 = d2 / (n - 1.0);    /* chisq per degree of freedom */
  252.  
  253.     *cov_11 = s2 * 1.0 / (n * (m_x * m_x + m_dx2));
  254.  
  255.     *sumsq = d2;
  256.   }
  257.  
  258.   return GSL_SUCCESS;
  259. }
  260.  
  261.  
  262. int
  263. gsl_fit_wmul (const double *x, const size_t xstride,
  264.           const double *w, const size_t wstride,
  265.           const double *y, const size_t ystride,
  266.           const size_t n, 
  267.               double *c1, double *cov_11, double *chisq)
  268. {
  269.  
  270.   /* compute the weighted means and weighted deviations from the means */
  271.  
  272.   /* wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i) */
  273.  
  274.   double W = 0, wm_x = 0, wm_y = 0, wm_dx2 = 0, wm_dxdy = 0;
  275.  
  276.   size_t i;
  277.  
  278.   for (i = 0; i < n; i++)
  279.     {
  280.       const double wi = w[i * wstride];
  281.  
  282.       if (wi > 0)
  283.     {
  284.       W += wi;
  285.       wm_x += (x[i * xstride] - wm_x) * (wi / W);
  286.       wm_y += (y[i * ystride] - wm_y) * (wi / W);
  287.     }
  288.     }
  289.  
  290.   W = 0;            /* reset the total weight */
  291.  
  292.   for (i = 0; i < n; i++)
  293.     {
  294.       const double wi = w[i * wstride];
  295.  
  296.       if (wi > 0)
  297.     {
  298.       const double dx = x[i * xstride] - wm_x;
  299.       const double dy = y[i * ystride] - wm_y;
  300.  
  301.       W += wi;
  302.       wm_dx2 += (dx * dx - wm_dx2) * (wi / W);
  303.       wm_dxdy += (dx * dy - wm_dxdy) * (wi / W);
  304.     }
  305.     }
  306.  
  307.   /* In terms of y = b x */
  308.  
  309.   {
  310.     double d2 = 0;
  311.     double b = (wm_x * wm_y + wm_dxdy) / (wm_x * wm_x + wm_dx2);
  312.  
  313.     *c1 = b;
  314.  
  315.     *cov_11 = 1 / (W * (wm_x * wm_x + wm_dx2));
  316.  
  317.     /* Compute chi^2 = \sum w_i (y_i - b * x_i)^2 */
  318.  
  319.     for (i = 0; i < n; i++)
  320.       {
  321.     const double wi = w[i * wstride];
  322.  
  323.     if (wi > 0)
  324.       {
  325.         const double dx = x[i * xstride] - wm_x;
  326.         const double dy = y[i * ystride] - wm_y;
  327.         const double d = (wm_y - b * wm_x) + (dy - b * dx);
  328.         d2 += wi * d * d;
  329.       }
  330.       }
  331.  
  332.     *chisq = d2;
  333.   }
  334.  
  335.   return GSL_SUCCESS;
  336. }
  337.  
  338. int
  339. gsl_fit_mul_est (const double x, 
  340.                  const double c1, const double c11, 
  341.                  double *y, double *y_err)
  342. {
  343.   *y = c1 * x;
  344.   *y_err = sqrt (c11) * fabs (x);
  345.   return GSL_SUCCESS;
  346. }
  347.